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Abstract. The use of spectral projection based methods for simulation of a stochastic system with discon- 
tinuous solution exhibits the Gibbs phenomenon, -which is characterized by oscillations near dis- 
continuities. This paper investigates a dynamic bi-orthogonality based approach with appropriate 
post-processing for mitigating the effects of the Gibbs phenomenon. The proposed approach uses 
spectral decomposition of the spatial and stochastic fields in appropriate orthogonal bases, while 
the dynamic orthogonality condition is used to derive the resultant closed form evolution equations. 
The orthogonal decomposition of the spatial field is exploited to propose a Gegenbauer reprojection 
based post-processing approach, where the orthogonal bases in spatial dimension are reprojected 
on the Gegenbauer polynomials in the domain of analyticity. The resultant spectral expansion in 
Gegenbauer series is shown to mitigate the Gibbs phenomenon. Efficacy of the proposed method is 
demonstrated for simulation of a one-dimensional stochastic Burgers equation with uncertain initial 
condition. 
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1. Introduction. Continuous advancements in digital technologies have resulted in ubiqui- 
tous use of computer simulators for investigation of many large-scale systems. The simulator 
predictions are often uncertain due to unkno-wn/poorly kno-wn physics, model parameters, 
initial and boundary conditions. Need and importance of uncertainty quantification in simu- 
lation predictions has already been emphasized by researchers in varied fields [34, 28, 29, 33]. 
Subsequently, research community have invested significant efforts to develop various uncer- 
tainty quantification methodologies (see [45] and references therein.). In essence, uncertainty 
quantification is an estimate of joint probability distribution of system responses conditional 
on probabilistic specification of uncertainties in the simulator model, parameters, initial and 
boundary conditions etc. 

1.1. Background. For further discussion, consider a probability space {Vl,J-,V) -where 
a; G is a set of elementary events, T is associated cr-algebra and P is a probability measure 
defined over J^. A generic stochastic partial differential equation (SPDE) on {Q,J^,V) is given 
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by 

— = C[u{x,t;uj);u}] (1.1) 

where u{x, t; uj) is a random variable, t G T is time, x ^ X C TZ'^ is a spatial dimension and C 
is an arbitrary differential operator of the form 

C[u{x,t;uj);uj] = — , (1.2) 

f{u{x, t; w); w) being flux. The initial condition of the system is specified using a random field 
u{x,0;uj), while, the boundary conditions are given by 

^(/3, t; Lo) = h{p, t; w); f^edX, w G (1.3) 

where i3 is a linear differential operator. Without loss of generality, the discussion in this paper 
is presented assuming u{x, t; u) to be a scalar function defined over a one-dimensional space, 
however, the discussion can be extended to a more generic case in straightforward manner. 

Owing to generality and ease of implementation, Monte Carlo methods are widely used for 
solution of (1.1) [37]. However, the Monte Carlo method typically requires a large number of 
samples for acceptable accuracy [37, 5]. For simulation of large scale systems, computational 
cost may render application of Monte Carlo methods for uncertainty quantification intractable. 
Stochastic spectral projection (SSP) based methods can provide a computationally efficient 
alternative to Monte Carlo methods with comparable accuracy. Under a generic condition 
of square integrability, u{x, t; co) can be represented as a parametric field in a Hilbert space 
T-L{T X X X [1]. The SSP methods use the tensor product representation of T-L^T x X x 
in appropriate Hilbert spaces Tii and with possible choices [44] 

n{Txx xF) = ni{TxF)xn2{x) (1.4) 

= ni{T)x'H2{X xF) (1.5) 

= ni{T xx)xn2{:F). (1.6) 

A Karhunnen-Loeve expansion [12] uses the tensor product (1.4), where u{x,t;u}) is spectrally 
represented in terms of orthogonal basis in 7i2{X) with temporally varying expansion coeffi- 
cients in TiiiT X F). A generalized Polynomial Chaos (gPC) expansion [50] is an example of 
tensor product (1.6) with orthogonal basis in 'H2{F) and the associated expansion coefficients 
in 'Hi{T X X). This paper investigates a bi-orthogonal method, which uses the gPC basis in 
H.2{^): while, dynamically orthogonal eigenfunction basis is used in 'Hi{T x X). 

One of the earliest exposition of a SSP-based method is the Homogeneous Chaos theory 
introduced by Wiener [47, 48], where random variables are expanded in terms of orthogo- 
nal Hermite polynomials. The expansion converges in mean square sense for any nonlinear 

functional [2]. Meecham and Jeng [27] used Homogeneous Chaos for turbulent mod- 
eling, however it was found that the convergence of chaos expansion is slow for turbulent 
field [38, 6]. Ghanem and Spanos [11, 10, 12] have proposed Polynomial Chaos (PC) method 
based on Homogeneous Chaos theory for uncertainty propagation in structural mechanics 
simulations. Subsequently, the method has been widely used for uncertainty propagation in 
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fluid flow simulations [21, 24, 26, 32]. Present state of the art for SSP methods is based on 
the Galerkin projection based generalized Polynomial Chaos method introduced by Xiu and 
Karniadakis [49, 50], which uses the representation of form (1.4). Using separability, 'Hi{J-) 
is decomposed in orthogonal polynomial chaos basis ^/^^(a;), where tpk{^^) belong to the Askey 
scheme of orthogonal polynomials [19] with weight function given by probability distribu- 
tion function of an uncertain input. Thus, using gPC, the dynamical response u{x, t; co) is 
spectrally represented as 

p 

u{x, t;uj) = Ui{x, t)'ipi{uj), (1.7) 
1=1 

where Ui{x,t) are dynamical expansion coefficients. See Mathelin et al. [25] and Najm [30] for 
extensive review of the research work in this field. 

Remark 1. Throughout this paper, the operator ' =' is used interchangeably to represent 
'equal to ' and 'approximately equal to '. Note that all the instances of ' =' in spectral repre- 
sentation defines the 'approximately equal to' operator. 

1.2. Motivation. Research work presented in this paper is motivated by the need to 
quantify uncertainty in simulation of systems with discontinuous solutions. In particular, 
the paper deals with stochastic systems defined using the hyperbolic conservation laws that 
invariably results in discontinuous solution irrespective of the initial conditions. Simulation 
of such system have attracted significant interest by research community and literature is 
rich with methods for resolution of discontinuities in deterministic simulations (see [22] and 
references therein for review of the methods). 

The use of spectral methods to hyperbolic SPDEs is known to pose following significant 
numerical challanges [31] 

1. numerical solution using finite difference schemes lead to spurious oscillations in ex- 
pansion coefficients and/or orthogonal basis 

2. spectral expansion of the solution of hyperbolic SPDEs lead to Gibbs phenomenon 
as discontinuities develop in solution [31]. The Gibbs phenomenon is characterized 
by [17, 14] a) slow convergence of spectral approximation of the solution at points 
away from the discontinuity; and b) 0(1) oscillations are observed in the solution near 
discontinuities that do not decrease with increasing N. 

There are recent research efforts that focus on addressing the issue of Gibbs phenomenon 
in gPC settings [31]. Wan and Karniadakis [46] have proposed a multi-element gPC (ME- 
gPC) method to mitigate the effect of Gibbs phenomenon (also see Lin et al. [23]). The 
ME-gPC method uses domain decomposition of the stochastic space with locally defined gPC 
basis, while the forward problem is solved for each subdomain. However, the computational 
cost of the ME-gPC method rises quickly as the number of subdomains increases. Poette et 
al. [36] have proposed an entropy based intrusive polynomial moment method for uncertainty 
propagation in non- linear systems with discontinuous solutions. The method reformulates the 
SPDE in terms of appropriately selected entropic variables that are bijection of the uncertain 
parameters, while, the gPC expansion of the entropic variables is defined through the Galerkin 
projection of the bijection. The method is implemented in two steps: in the first step, the 
resultant stochastic Galerkin system is discretized using a finite volume approach, and up- 
winded Row scheme is used for numerical solution. In the second step, gPC expansion of the 
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entropic variables is calculated so that the entropy of the system is minimized. The method is 
demonstrated for solution of the stochastic Burgers equation. The method is computationally 
efficient than the ME-gPC, however, the requirement for minimization of the entropy at each 
time step results in the computational cost higher than the gPC method. Tryoen et al. [43] 
uses Roe-type solver for intrusive Galerkin projection of uncertain hyperbolic systems. The 
solver uses Dubois and Mehlman entropy corrector that avoid certain entropy violating solu- 
tions. Sargsyan et al. [40] have proposed a two step method for uncertainty quantification of 
systems with discontinuous response. In the first step, a limited number of simulation runs 
are used to infer the shock location using the Bayesian inference. In the second step, localized 
gPC basis is defined in the region of analiticity, while intrusive Galerkin projection is used 
to propagate the uncertainty. In an alternate non gPC setting, Chantrasmi et al. [4] have 
proposed a Pade-Legendre interpolant based approach, where, simulation runs are obtained at 
predefined Gauss-Legendre quadrature nodes and the resultant discontinuous system response 
is reconstructed using the Pade interpolation. Despite of these recent progresses, resolution 
of Gibbs phenomenon in application of SSP methods to stochastic hyperbolic systems is an 
area of current research interest [40]. 

1.3. Proposed Method. Existence and resolution of Gibbs phenomenon for determinis- 
tic functions is extensively investigated in the literature [17, 14, 16]. Spectral expansion of 
a function using global orthogonal basis is contaminated by the presence of discontinuities. 
However, expansion coefficients contain enough information to recover the function with high 
accuracy using appropriate post-processing. Gibbs phenomenon can be completely resolved by 
reprojecting the partial sum on Gibbs complementary basis [14]. Gottlieb et al. [15, 13, 14, 16] 
have shown that exponentially convergent approximation can be obtained at point values in 
subinterval of analyticity by reprojecting the partial sum on orthogonal Gegenbauer polyno- 
mial basis. The post-processing method can be extended to stochastic functions, provided, 
the spectral expansion is available in terms of orthogonal basis in spatial dimensions. 

Note that the gPC expansion (1.7) uses orthogonal basis in (-^) while corresponding 
coefficients Ui{x,t) are functions in 'Hi{T x X), though, the discontinuity reside in spatial 
dimension. Thus, enough information is not available in Ui{x, t) to resolve Gibbs phenomenon 
using Gegenbaeur reprojection method. Resolution of Gibbs phenomenon proposed in this 
paper is based on the observation that exponential accuracy can be recovered if u{x, t; uj) is 
expanded in terms of orthogonal bases in spatial dimension, and subsequently reprojected on 
an approximate Gibbs complementary basis. 

Tagade and Choi [42, 41] have proposed a dynamic bi-orthogonal field equations (DBFE) 
method for solution of SPDEs in the context of Bayesian calibration. This paper proposes a 
method for solution of hyperbolic SPDEs using DBFE, that exploits orthogonal decomposition 
of the spatial dimension to develop a post-processing approach for mitigation of the Gibbs 
phenomenon. Consider a generic Karhunnen-Loeve expansion of u{x, t; u) 



i=l 

where u{x,t) is the mean, Ui{x,t) are eigenfunctions which form complete orthogonal basis 
in T-Li{X) at a given time step t, while Yi{t;uj) are independent zero-mean random variables. 
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Using polynomial chaos expansion of Yi{t]uj) in (1.8), bi-orthogonal expansion of u{x,t]uj) is 
given by 

N P 

u{x,t;Lj) =u(2;,i) + ^^y;(t)V'p(w)Mi(x,t). (1.9) 
1=1 p=i 

A Dynamic Orthogonality (DO) condition [39] is used to derive closed form solution of evolu- 
tion equations for u{x, t), Ui{x, t) and Yp{t). However, the resultant approximation oiu{x, t; oj) 
in (1.9) using the solution olu{x,t), Ui{x,t) and Yp(t) exhibits the Gibbs phenomenon. A 
Gegenbauer reprojection based post-processing method is proposed to mitigate the effect of 
Gibbs phenomenon by reprojecting the eigenfunctions Ui{x,t) on the Gegenbauer basis. The 
proposed method is demonstrated for solution of a stochastic one-dimensional Burgers equa- 
tion with uncertain initial conditions. 

The rest of the paper is organized as follows. In section 2, mathematical formulation 
of the method is provided. Proposed post-processing method for Resolution of the Gibbs 
phenomenon using dynamic bi-orthogonality based method is discussed in section 3. Section 
4 provides numerical results for Burgers equation to demonstrate efficacy of the proposed 
method. Finally, paper is summarized and concluded in section 5. 

2. Dynamically Bi-orthogonal Field Equations (DBFE). Definition of the following op- 
erators will be used in the remainder of the paper. 

Definition 1. Inner product on 'Hi{T x X) for a given t is defined as 



{u{x,t;uj),v{x,t;uj)) = / u{x,t;uj)v{x,t;u})dx. (2.1) 

Jx 

Similarly, inner product on 7^2 (-^) is defined as 

{u{x,t;u)),v{x,t;uj))^ = / u{x,t;u))v{x,t;uj)dV{uj). (2.2) 



The expectation operator on T-L2{J^) is defined as 

u{x,t) = E'^[u{x,t;uj)]= / u{x,t;uj)dV{oj). (2.3) 

Further using definition of the expectation operator, the covariance is defined as 

Cu,v = E'^[u{x, t; uj)v{x, t; uj)]. (2.4) 



The DBFE method uses the bi-orthogonal expansion (1.9) of u{x,t;uj) in (1.1). How- 
ever, independent evolution equations for the mean u{x, t), the eigenfunctions Ui{x, t) and the 
corresponding coefficients Yp{t) can not be obtained concurrently using the SPDE in (1.1). 
Well posed evolution equations for the unknown quantities can de derived by imposing the 
additional dynamic orthogonality (DO) condition [39]. The DO condition is specified as [39] 

(^u.ix,t),^^^^^^ = 0, Vz, j = l,...,7V, (2.5) 
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where N is the number of eigenfunctions retained in the expansion (1.9). Note that the DO 
condition constrains the dynamic evolution of the eigenfunction basis such that the orthonor- 
mality of the eigenfield is retained at all time steps [39]. 

Using the DO condition and bi-orthogonal expansion of u{x,t;u}), SPDE (1.1) can be 
reformulated into a set of + 1 PDEs and N x P ODEs as follows. 

Proposition 1. Using the DO condition, dynamic evolution equations ofu(x,t), Ui{x,t) and 
Yp{t) are given by 



du{x, t) 

dt 
dui{x, t) 
di 



E'^ [C[u{x,t,;uj);uj]] (2.6) 

N 

Cy},^{E^ [C[u{x, t, oo)-iv]Yj (t; iv)] 



N 



J2 {E'^iCHx, t, u;y,oj]Yj{t; u;)],Uk{x, t))^ } (2.7) 



fc=i 



dY;{t) _ 1 



{{C[u{x,t;uj);uj] - E'^[C[u{x,t;uj);uj]],Ui{x,t)) ,ipp{uj))^ . {2.i 



Proof. Proof of (2.6) and (2.7) is given by Sapsis and Lermusiaux [39], which is briefly pre- 
sented here for completeness, while, (2.8) is derived here by introducing the bi-orthogonality. 
Proof of (2.6) and (2.7). Use a generic KL expansion (1.8) in (1.1) to obtain 

i=l i=l 

Application of the expectation operator to (2.9) gives (2.6). 

Multiply (2.9) by Yj(t;ijj) and apply the expectation operator to obtain 

^ ^ av (r t) 

E^iM!)y^,(,)^^(^>*) + J2Cy.mAt)^^ = [u{x,t;ujy,u] Y,{t;u)] . (2.10) 

i=l i=l 

Multiply (2.10) by Uk{x,t), take inner product and apply the DO condition to get 

CdY,(t)^..^ = {E'^[C[u{x,t;u;y,uj]Yjit;u;)],Uk{x,t))^. (2.11) 

Use (2.11) in (2.10) to obtain 

ECy.wy.w^^^^ = ^"[^K^>i>^);M^.-(*;^)] (2-i2) 



4 = 1 



N 



^(^"[£[n(x,t,u;);a;]yj(t;a;)] ,Uk{x,t))^Ukix,t), (2.13) 



k=l 
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which can be written in the matrix as 

U = E-^D, (2.14) 

where S is the covariance matrix with ij^^ element Sjj = Cy.y^. . 

Multiply both sides of (2.9) by Uj{x,t), take inner product and use the DO condition and 
orthonormality of eigenfunctions to get 

/ du(x,t) , ,\ dYi(t:u}) / N 1 / N> 

y ,Ujix,t)J^+ ^ = {C[u{x,t;u;y,Lo],Uj{x,t))^, (2.15) 

which on application of expectation operator gives 



(^^^^^^,u^{x,t)^ =E'^ [{£[u{x,t;ujy,u],Uj{x,t))^] . 



(2.16) 



Using (2.16) in (2.15) gives 

dYi{t;Lo) ^ 
dt 



([£ [{x, t; a;); - E'^ [C [u(x, t; a;); w]]] , mix, t)) ^ . (2.17) 



Proof of (2.8). Basic conditions of KL expansion ensures that Yi{t]uj) are square inte- 
grable random variables, thus, according to Cameron and Martin theorem [2], Yi{t;oj) can be 
approximated to arbitrary accuracy using spectral expansion in terms of orthogonal basis in 
T-l2{J^)- Hence, Yi{t]uj) can be spectrally represented using P terms as [50] 

p 

Yi{t;u) = Y,Y;mp{co), (2.18) 
p=i 

where ipp{ijj) are orthogonal basis in {^)- Differentiating (2.18) with respect to t gives 

dY.jt-u:) ^dY;{t) 

-^^ = l^^^M^l (2.19) 

p=i 

which on Galerkin projection provide 

dYi^t) \—k^^M^)) 



n 



Having (2.19) and (2.20) in (2.17) results in (2.8). 



(2.20) 
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2.1. Initial and Boundary Conditions for DBFE. The initial conditions for DBFE are 
defined by specifying the mean u(x,0), the eigenfunctions Ui{x,0) and the coefficients Yp{0), 
which are given by considering the bi-orthogonal expansion of the initial condition of SPDE 
(1.1) as 

N P 

«(x,0;w) =n(x,0) + y; (OK (x, 0)^-^(0;). (2.21) 

1=1 p=i 

Applying the expectation operator, initial condition for the mean field is given by 

u{x,0) = E'^ [u{x,0;u)]. (2.22) 

The initial condition for eigenfunctions, Ui{x,0), is given by solution of the Fredholms' 
integral equation of the second kind [3] 

Cuixt,o),ui^2,o)M^i^^)dxi = X^Ui{x2,0) (2.23) 



X 



where Cy_(^xi,o),u{x2,o) is the covariance function of u{x,0;uj) and Aj are the eigenvalues. A 
Galerkin projection based method is used to solve (2.23) numerically [18]. For a Gaussian 
process, all the coefficients Yp{0) are zero except 



>1(0) = VAi, (2.24) 

thus defining Yi{0;uj) as normal variables with variance Aj. For a generic non-Gaussian case, 
the initial expansion coefficients Yp{0) are given by 



{u{x,0;uj) - u{x,0),Ui{x,0));^ ,'iljp{uj) 

Boundary conditions for DBFE are given by using the generic KL expansion of h{l3,t;uj), 
which specifies the boundary conditions for the SPDE in (1.1), as 

N 

h{P,t;uj) = h{P,t) + J2y^{t■,0J)ui{P,t). (2.26) 

i=l 

Applying the expectation operator on (2.26), boundary condition for the mean field is given 
by 

B[u{P,t;uj)]^^g^ = h{^,t). (2.27) 
Multiply (2.27) by Yj{t;uj) and apply the expectation operator 

N 

[h{fi,t;uj)Yj{t-u:)\ = Y,CY,{t)Y,{tMP,t), (2.28) 

i=l 

which can be specified in a matrix form as 

u = S^^E, (2.29) 

where S is the covariance matrix. 
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3. Resolution of the Gibbs Phenomena using Dynamically Bi-orthogonal Field Equa- 
tions . Key challenges involved in the numerical implementation of the DBFE for SPDEs 
with discontinuous solutions are: (a) numerical evaluation of the mean u{x, t) and the eigen- 
functions (using (2.6) and (2.7)) results in the development of spurious oscillations as dis- 
continuities evolve in the solution; (b) the resultant bi-orthogonal expansion (1.9) exhibits 
the Gibbs phenomenon characterized by the oscillations near the discontinuity location and 
the slow convergence away from the discontinuity location. This paper proposes a two step 
approach for the application of DBFE to the stochastic systems with discontinuous solution, 
that exploits the orthogonal decomposition of the spatial field in the DBFE approach. In the 
first step, extension of the existing schemes is proposed to derive the non-oscillatory numerical 
scheme for DBFE. In the second step, a Gegenbauer reprojection based method is proposed 
to mitigate the effects of the Gibbs phenomenon. In this section, the proposed approach is 
described in detail. 

3.1. Numerical Scheme for Non-oscillatory Solution. Due to discontinuities, numerical 
solution of (2.6)-(2.8) is expected to contain spurious oscillations with reduced accuracy in 
u{x,t) and Ui{x,t). Total Variation Diminishing (TVD) finite difference schemes are widely 
used in the literature to obtain non-oscillatory accurate weak solutions to hyperbolic partial 
differential equations [22]. This paper demonstrates extension of a TVD scheme for numerical 
solution of (2.6)-(2.8), however, note that other numerical schemes can similarly be extended 
for DBFE. 

For notational convenience, this subsection defines u^j{uj) as an approximate value of u{x = 
Xj,t;uj) at the grid point Xj = jAx and present time t. Consider a numerical scheme applied 
to SPDE (1.1) for a given u as 

duUuj) 1 / - - \ 

' - (3.1) 



dt Ax ^^+2^ ^ 2 

where /j_,_i(t<j) is an interpolated flux inside a cell [rrj, Xj+i]. A {2k + 1) point numerical 
scheme is defined by using k points on either side of Xj for interpolation, thus, 

f^^iu) = /(4_,+i(a;),...,4+,(a;)), (3.2) 

such that [22] 

f{u,...,u) = f{u). (3.3) 
fi-li'^) is also defined in a similar manner. The TVD scheme (3.1) can be extended to DBFE 

J 2 

by using 

C[u{x,t;ujy,u^] = - (3-4) 

for numerical solution of (2.6)-(2.8). The resultant numerical scheme retain non-oscillatory 
property in u{x,t) and Ui{x,t). 

In the present paper, the DBFE method is implemented using the first order central 
scheme proposed by Kurganov and Tadmor [20] (central KT scheme), for which numerical 
flux is defined as 

(w) = -I [/(^+i(u;)) + /(^(^))] + I [a,.+ i (^) (n*+i(^) - n*(a;))l , (3.5) 



2 
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where a-,i(oj) is the maximum local speed at x,-. The gPC expansion of a-,i(uj) is given by 

aj+i(w) = 'V'p(^), (3.6) 



where ctp ^ are gPC expansion coefficients. Use (3.6) in numerical flux (3.5) to obtain 



TV 



j=l 



1 



--{i^- [/(^+i(a;))n(t;a;)] [/(n* (a;))n(t; ^)] } 

+i{(Mf'^+^[n(x,+i,t) 



N 



1 



f.M^) - E- ! (c) = --{/(^^,(a;)) - i?- + /(n*(a;)) - i^^ [/(^(u;))] } 



Af P ^ 



(3.7) 



where ^jj+i = '5^j(i; '^)«j+i (w) and ^. ^^ = yj(t; (ii;)yfc(t; (ii;)a^._^i (c^;). 

The numerical flux (3.7) is used in (3.1) to obtain the first order central difference scheme 
for numerical solution of (2.6)-(2.8). 4*^^ order Runge-Kutta method is used for time inte- 
gration in (2.6) and (2.7), while Euler method is used to solve ODE (2.8). Gauss-Legendre 
quadrature is used to calculate inner products in spatial dimensions, while, Monte Carlo 
integration is used for inner products in stochastic dimensions. 

3.2. Post-processing using Gegenbauer Polynomial. Let the solution of the SPDE in 
(1.1) is obtained using the DBFE method till time t = T such that the discontinuities are 
developed. For a given uj, let [a{u)), b{uj)] denote an interval of analyticity for u{x, T; u), where 
a random variable C{x; co) can be defined such that —1 < C{x; uj) < 1. Gibbs phenomenon can 
be resolved in the interval [a{uj), b{u;)] by reprojecting the partial sum (1.9) on a suitable Gibbs 
complementary polynomial basis [16]. In the present paper, the partial sum is reprojected on 
a Gegenbauer Polynomial, which is defined for a C{x; u) as follows. 

Definition 2. [16] The Gegenbauer polynomial C^{C,{x]Oj)) is a polynomial of degree n which 
is orthogonal over [—1, 1] with weight function (1 — ('^{x;uj))'^~2 for A > 0. The orthogonality 
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is given by 



where 



(1 - Cix;u))^-^C^:{Cix;u;))Ci{ax;io))dC = K6m,. 



V^c„^(i)r(A + i) 



with 



r(A)(m + A) ' 
r(n + 2A) 



m!r(2A) 



C^^{({x;u})) can be estimated using a recurrence relationship 



(3.8) 
(3.9) 
(3.10) 



(n + l)C^+i(C(x; u)) = 2(A + n)C(x; u;)C7^(C(x; a;)) - (2A + n - l)C^i(C(x; u)). (3.11) 



Using orthogonality of Gegenbauer polynomial, u{x,T;uj) can be represented in terms of 
exponentially convergent Gegenbauer expansion series truncated at M terms as 



M 



uix 



(3.12) 



where ({x;u}) = ^^7^|y^; e('^) = ^"(i^) ^^^-^ _ &(aj)+a(a;) ^ Qoefficients diico) are given by 



(3.13) 



However, u{x{(),T;uj) is not available for calculation of ui{u}). Using bi-orthogonal expansion 
(1.9) in (3.13), approximate Gegenbauer expansion coefficients gi{uj) are given by 



91 i^) 



1 



1 - C{x- u)f~-^Ct{ax: uj))u{x, T)dC 



N P 1 

Y,Y.^;iT)i;,{u) / {l-e{x;u))'-2C^{ax;u))u,{x,T)dC 
i=i p=i -^-1 



(3.14) 



Note that (3.14) uses a bi-orthogonal expansion of u{x,T;lj), thus, gi{u}) is an approximation 
of ui{uj). gi{uj) can be expanded in gPC basis as 



where 



{gi{u),'4jp{u))^ 



(3.15) 



(3.16) 
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Coefficients are calculated using Gaussian quadrature by evaluating (3.14) at collocation 
points, while, the integrals in (3.14) are calculated using Gauss-Gegenbauer quadrature. Using 
g^p, u{x, T; to) is approximated as 

AI P 

u{x,T;uj) = ^^gl,iljp{oj)CHax;oj)). (3.17) 
1=1 p=i 

The partial sum (3.17) approximates u{x,T;uj) with exponential accuracy in the subinterval 
[a{Lo), b{ui)]. Note that determination of domain [a(w), b{uj)] requires location of discontinuity. 
An edge detection method [8, 9] can be used to locate points of discontinuities. 

4. Numerical Example: ID Burgers Equation. Efficacy of the proposed method is inves- 
tigated for solution of a stochastic one-dimensional Burgers equation with uncertain input con- 
ditions. Note that previous work has utilized deterministic Burgers equation extensively [51], 
while, there have been recent investigation on solution of stochastic Burgers equation [35, 36]. 
Coupled with the ease of implementation, stochastic Burgers equation provide an appropri- 
ate context to investigate the efficacy of the proposed method for mitigation of the Gibbs 
phenomenon. 

A one dimensional inviscid stochastic Burgers equation is given by 

where x G [—1,1]. In the present paper, the stochastic Burgers equation is investigated for 
uncertain initial condition given by 

u{x,0;uj) = sG{x;u}), (4.2) 

where G{x;uj) is a Gaussian process with mean 

u{x, 0) = Ub — tan"""^ (x — xq) , (4.3) 

where Ub and xq are user defined constants, and the covariance function 

Cu{xuO)M^2,o) = exp (-A | xi - 2:2 |) , (4.4) 

where a is standard deviation and A is the correlation length. In this paper results are 
presented for s = 0.1. The boundary condition is specified using 

u{x, t; oj) = u{(3, 0; w). (4.5) 

4.1. Monte Carlo Method for Stochastic One-dimensional Burgers Equation. To in- 
vestigate accuracy of the proposed DBFE method, its numerical results are compared with 
the Monte Carlo simulation. The Monte Carlo method for a stochastic one-dimensional Burg- 
ers equation is applied by solving the deterministic Burgers equation with initial condition 
defined using samples of u{x,0;uj). In this section, results are presented for a test case with 
Ub = 0,xq = 0. For each sample, explicit fourth order Runge-Kutta scheme is used for time 
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evolution while the first order central KT scheme [20] (see equation (3.5)) is used in spatial 
dimension. The numerical solutions are obtained for Ax = 0.01 and At = 10"^. Total 10000 
samples are collected for the Monte Carlo simulation. Figure 4.1(a) shows several realizations 
of initial conditions, where the mean is shown with emphasis, while respective solutions at 
time t = 1.1s are shown in Figure 4.1 (b). Shock formation is observed for all the samples, 
however, shock location varies depending on the initial condition. 




(a) 




(b) 



Figure 4.1. Solution of SPDE (1.1) using Monte Carlo method. Figure a) shows samples of initial condi- 
tion; figure h) shows solution of samples corresponding to the initial conditions in a). 



4.2. Application Results. Use the bi-orthogonal expansion of u{x,t;uj) in (4.1) to obtain 

2 

1=1 p=l 



dx 2 

N N P P 



E E E E Y;{t)Y^{t)i;,{u)^l.,{u:) 

i=l j = l p=l q=l 



d f Ui{x, t)uj{x, t) 
dx 



■ (4.6) 



The differential operator (4.6) is used in (3.7) to derive the DBFE governing equations (2.6- 
2.8). 

Initial condition for the mean u(x, t) is specified using (4.3), while the initial condition for 
Ui{x,t) and Yp{t) are given by solution of the eigenvalue problem for the covariance function 
(4.4). Note that the eigenvalues for the covariance function (4.4) decreases rapidly with only 
first three eigenmodes dominant, thus, = 3 suffices for the accurate approximation of the 
uncertain initial condition. However, to account for a highly nonlinear nature of the Burgers 
equation, higher value of N may be required. Figure 4.2 shows eigenvalues and first four 
eigenfunctions of the uncertain initial condition with a = 0.5 and A = 1.0. 
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Figure 4.2. Eigenfield for uncertain initial condition, a) eigenvalues and h) first four eigenfunctions 



4.2.1. Non-oscillatory Solution using the Proposed Numerical Scheme. Efficacy of the 
numerical scheme described in the section 3.1 to obtain non-oscillatory solutions for mean and 
eigenfunctions is investigated by implementing the DBFE method using the first order central 



KT scheme (3.7), central difference scheme (terms involving a^^'^ , Mp''^'^ and Tj 

neglected in (3.7)) and the first order central KT scheme in mean (terms involving aj"^^ and 

M-^' ^ are retained only for the first term in (3.7)). 

Figure 4.3 shows comparison of results for (a) mean and (b) first eigenfunction for = 3. 
With central difference scheme, oscillations are observed in both mean and eigenfield. When 
central KT scheme is used only in mean field and central difference for eigenfield, oscillations 
in mean are comparatively reduced. However oscillations are still present in eigenfield, as 
can be observed for first eigenfunction. When full central KT scheme is used, oscillations in 
both mean and eigenfield are resolved. Note that only averaged statistics in terms of expected 
values and covariance functions affect time evolution of mean and eigenfields, hence, no special 
treatment is necessary for random field. 

4.2.2. Characteristics of DBFE Solutions. The characteristics of the numerical solution 
for DBFE before post-processing is investigated. In order to investigate effect of truncation 
of KL expansion. Burgers equation is solved for = 3, = 5 and N = 7. To compare 
results with Monte Carlo simulations, u{x,t;u}) is reconstructed for every sample using (1.9). 
Accuracy of the DBFE method is compared with Monte Carlo using relative error in Li{il.) 
norm, which is given by 



are 



(x,t) 



!n!&x u{x,t,uj)- [ui,{x,t)+Y,i=iT.p=iY^{t)ui{x,t)'ipp{u)\ dxdV{uj 



fnlsx I u{x,t,uj) I dxdV{uj) 



(4.7) 
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Figure 4.3. Figure shows comparison of numerical solution for a) mean and b) first eigenfunction, obtained 
by implementing the proposed DBFE method using the central difference scheme, first order KT scheme m mean 
and full first order central KT scheme. The results are obtained using the first three eigenmodes (N=3) m bi- 
orthogonal expansion. 



Remark l.Note that the Li-error in (J^.l) is defined for a small domain 5x, and thus is not 
a point error. In the present paper, the Li-error at a point is defined using the adjoining two 
cells. 

Figure 4.4 shows the Li-error of the DBFE method relative to the Monte Carlo method. 
Li(u}) error of the order of 10~^ is observed away from the shock location, however, near the 
shock location, error of the order of 10~^ is observed that does not reduce significantly with 
increasing N. 

Figure 4.5 (a) shows the first eigenfunction at t = 1.1s for = 3, 5, 7. The eigenfunction 
shows periodic oscillations near shock, with number of modes increasing with N. Figure 
4.5 (b) shows evolution of variance of different modes over time for N = 7. Note that the 
modes with negligible variance at t = become dominant over time evolution. This effect can 
be significantly observed for the 5^^ eigenmode, which has negligible variance at t = but 
becomes the second most dominant mode over the time evolution. Thus, it is necessary to 
use higher number of modes, even though the modes show negligible variance at i = 0. 

Figure 4.6 shows comparison of the 90% confidence bound for the Monte Carlo and the 
DBEE method. Note that the confidence bound shows significant disagreement near shock 
location for all the test cases. The discrepancy in the confidence bound is due to the oscillations 
in the samples near the shock location, characterizing the Gibbs phenomenon for each sample. 

The key observations from the results presented here can be summarized as: a) numerical 
solution of the DBFE exhibits oscillations near the discontinuity, b) Li-error is maximum 
at the shock location, c) The amplitude of the maximum Li-error does not decrease by in- 
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Figure 4.4. Figure shows error for DBFE method. The error is calculated against the Monte Carlo 

method using 10000 samples. 




creasing the number of eigenfunctions, N . These observations are characteristics of the Gibbs 
phenomenon, indicating the need for appropriate post-processing. Efficacy of the proposed 
post-processing to mitigate the effect of the Gibbs phenomenon is demonstrated in the fol- 
lowing. 
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Figure 4.6. Figure shows comparison of 90% confidence hound. 

4.2.3. Mitigation of Gibbs Phenomenon using Post-processing. The Gegenbauer repro- 
jection based approach, outhned in section 3.2, is used for post-processing to mitigate the effect 
of the Gibbs phenomenon. The bi-orthogonal expansion of u{x, t; ui) is reprojected on 7^^ order 
Gegenbauer polynomials, while total 7 coefficients are used in the expansion (3.17). Integrals 
involving the Gegenbauer polynomials are numerically evaluated using Gauss-Gegenbauer 
quadrature with 100 nodes, whereas, the integrals involving the polynomial chaos basis are 
numerically solved using Monte Carlo integration. In the present paper, shock is assumed to 
be located at a point with highest slope where u = line is crossed. However, the proposed 
method can be implemented using any edge detection method for shock localization. 

Resultant samples after post-processing using reprojection method are compared with 
samples without post-processing in Figure 4.7 (the first six samples are shown in the figure). 
The respective Monte Carlo samples are also shown in the figure. Resolution of the Gibbs 
phenomenon can be observed for individual samples. For all the samples, oscillations are 
observed for approximations of u{x, t; oj) obtained using DBFE method. Magnitude of oscil- 
lations is low for samples near mean and increases away from mean. The oscillations show 
three distinct behaviors: (1) when shock is located at x > 0, oscillations are observed on left 
side of shock location; (2) when shock is located at x < 0, oscillations are observed on right 
side of the shock location; and (3) when shock is located near x = 0, oscillations are observed 
on both the sides with low magnitude. 

Figure 4.8 a) shows comparison of 90% confidence bound of DBFE method after post- 
processing with Monte Carlo simulation. Using post-processing, agreement with Monte Carlo 
simulation is significantly improved. Figure 4.8 b) shows Li error after post-processing. Com- 
paring with Figure 4.4, it can be observed that the error is reduced at all locations by the 
order of lO'^. At mean shock location, Li error is of the order of 10 . 

Figure 4.9 shows comparison of first four moments of DBFE solution with the Monte 
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Figure 4.7. Figure shows comparison of DBFE method with Monte Carlo samples. Comparison without 
and with post-processing are shown for first six samples. 
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Figure 4.8. Figure shows a) confidence hounds using post-processing b) Li error 



Carlo simulation. The comparison is shown for both without and with post-processing cases. 
For the higher moments, agreement with the Monte Carlo method improves significantly 
after the post-processing. Mean of the solution obtained using the DBFE method without 
post-processing matches well with the Monte Carlo method (shown in Figure 4.9 (a)). Note 
that the mean is not contaminated by the Gibbs phenomenon, resulting in the good match 
between DBFE method without post-processing and the Monte Carlo method, while, the 
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subsequent post-processing retains the accuracy of the mean. Figure 4.9 (b) shows the variance 
as a function of spatial location x. The variance is low in the region away from the shock 
location, whereas, significant increase in the variance is observed at mean shock location. Good 
agreement between the Monte Carlo and the DBFE method is observed in the region away from 
the shock location, however, nontrivial difference is observed without post-processing near the 
mean shock location, with DBEF predicting higher variance than the Monte Carlo. The high 
variance in DBFE predictions results due to the oscillations near shock location characterizing 
the Gibbs phenomenon, which are mitigated using the post-processing, providing the close 
agreement with the Monte Carlo method. Figure 4.9 also shows comparison for (c) skewness 
and (d) kurtosis. Skewness is zero in a region away from the shock location, indicating 
symmetric probability distribution. In a region near shock, skewness is negative in the left 
side of the mean shock location and positive on the right, with significantly high absolute 
value near the mean shock location x = 0. Similar results are obtained for kurtosis (Figure 
4.9 (d)) where high values are obtained near shock. It can be observed that DBFE method 
has predicted the trend correctly for skewness and kurtosis, while the close agreement with 
the Monte Carlo method is obtained after post-processing. 

4.3. Computational Efficiency of the Proposed Method. The computational cost of 
the numerical implementation of the proposed DBFE method is compared with the Monte 
Carlo and the gPC method. See Pettersson et al. [35] for the governing equations of the 
gPC method applied to stochastic Burgers equation. The method is implemented with S*""^ 
order gPC basis. Figure 4.10 shows the 90% confidence bound for the gPC method, which 
is compared with the Monte Carlo method. Comparison for the mean is also shown in the 
figure. The mean obtained using the gPC method agrees well with the Monte Carlo method, 
however, the confidence bound shows oscillation near the mean shock location, demonstrating 
the existence of the Gibbs phenomenon for solution of the SPDE (1.1) using the gPC method. 
Note that the similar observations are reported earlier in the literature for solution of the 
stochastic Burgers equations using the gPC method [35] . 

To investigate the computational efficiency of the proposed method, the Monte Carlo 
method, the gPC method and the DBFE are implemented on a desktop computer with Intel 
Core 15 CPU. Total 10000 samples are used for the Monte Carlo method, while, 3'^'^ order poly- 
nomial chaos basis is used for the DBFE and the gPC method. Time required for simulations 
is obtained using the FORTRAN intrinsic routine cpu-time. The comparison of CPU time 
is shown in Figure 4.11. Both the gPC and the DBFE methods provide considerable com- 
putational speed-up over the Monte Carlo method, with the computational time requirement 
increasing with the number of eigenfunctions used. Even for N = 7, speed up of the order 
of 3 is obtained as compared to Monte Carlo method. The computational cost of the DBFE 
method is lower than the gPC method, with the ratio of the CPU time for gPC and the DBFE 
decreasing with increasing A^. Although the post-processing increases the computational cost 
of the DBFE method, the DBFE method still remains computationally efficient with the CPU 
time approaching that of the gPC with A^. Note that the computational cost of the proposed 
method depends on the complexity of the differential operator C[u{x,t;u});uj], however, the 
computational cost of the post-processing remains constant irrespective of the complexity of 
the SPDE. 
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Figure 4.9. Figure shows comparison of first four moments for Monte Carlo and DBEF. a) mean b) 
variance c) skewness d) kurtosis 



4.4. Choice of Gegenbauer Parameters. The proposed reprojection method requires 
specification of parameter A of C^{x) and maximum number of expansion terms M. Total 
error in the Gegenbauer reprojection method eminates from error due to approximation of 
a function using a finite Gegenbauer series expansion and roundoff error in computation 
of Gegenbauer polynomials C^(x). Roundoff error in evaluation of C^{x) increases with 
increasing A and n [7], rendering very large A and M unsuitable. 

To investigate the effect of choice of parameters of Gegenbauer polynomial, proposed post- 
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Figure 4.10. Figure shows comparison of the confidence bound obtained using the gPC method with the 
Monte Carlo method. Comparison for mean is also shown in the figure. 
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Figure 4.11. Figure shows comparison of computational time for simulation. 10000 samples are used for 
the Monte Carlo method. JT"^ order Hermite polynomials are used as basis for DBFE and gPC 



processing method is applied using different values of A and M. Figure 4.12 shows minimum 
value of A for which Li error is greater than the tolerance for different N . For given M, 
limiting value of A decreases with A^. In the present study, TOL = 0.1 is used, note that 
since u{x, t; uj) is linear on the both side of shock, gi{uj) is dominant for all the samples, while 
values of higher coefficients are very low. The higher coefficients control non-linearity near 
shock. For M = 3, highest order of Gegenbauer polynomial is quadratic in nature, which 
can not properly capture non-linearity near shock. Since the values of coefficients decrease 
rapidly, roundoff errors dominates for higher Gegenbauer polynomials resulting in resulting 
in oscillations near shock discontinuity. For the test case presented in this paper, choice of 
A = 7.0 and M = 7 provide trade-off between requirement for capturing non-linearity near 
shock and lower round-off error. Although the choice of parameters is heuristic in nature, a 
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robust optimum parameter selection method can be developed based on already established 
techniques for deterministic solutions [7]. 
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Figure 4.12. Figure shows limiting value of X for given tolerance. 

4.5. Robustness of the Method. Robustness of the proposed method is investigated for 
following covariance functions: 1) squared exponential C{xi,X2) = cr'^exp [—{xi — X2)^), 2) 
triangular C{xi,X2) = (1.0 — t{\ xi — X2 \))exp{— \ xi — X2 |) and 3) uniformly modulated 
C{xi,X2) = a'^exp {—{xi + X2))exp{— \ xi — X2 |). Figure 4.13 shows Li error for different 
values of cr^ when squared exponential covariance function is used to represent uncertainty 
in initial conditions. Li error increases exponentially with increasing cj^, while, the error 
decreases after post-processing for all cj^. For high cj^, coefficients of higher modes of eigenfield 
become significant, necessitating the use of high A^. Thus, error due to truncation of spectral 
expansion at terms dominates the Li error for high cr^, resulting in lesser improvement 
using post-processing. 

Figure 4.14 shows comparison of 90% confidence bound obtained using MC and DBFE 
for squared exponential, triangular and uniformly modulated covariance function. Confidence 
bounds for DBFE method agrees well with MC, demonstrating robustness of the proposed 
method to resolve Gibbs phenomenon for different choices of covariance functions. 

5. Concluding Remarks. This paper have proposed a dynamic bi-orthogonality based 
spectral projection method for uncertainty quantification for systems with discontinuous solu- 
tions. The proposed approach is implemented in two steps: in the first step, input uncertainty 
is propagated to the system response using DBFE method, while, in the second step, effect of 
the Gibbs phenomenon is mitigated by reprojecting the mean and the eigenfield on the Gegen- 
bauer polynomials. Efficacy of the proposed method have been investigated for solution of a 
one-dimensional stochastic Burgers equation. The numerical results presented in this paper 
have demonstrated the ability of the proposed post-processing method to mitigate the effect 
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Figure 4.14. Figure shows comparison of confidence bound for different test cases. 

of the Gibbs phenomenon as discontinuities develop in the solution. Note that the proposed 
method does not require a-priory knowledge of the shock location, thus, a generic implemen- 
tation for a variety of applications can be achieved by extending any numerical scheme to the 
DBFE, as demonstrated in this paper. The DBFE method is found to be computationally 
efficient than the gPC method, thus, the method becomes an attractive alternative to the 
gPC for solution of SPDEs. 
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